
rm(list=ls())

library(tidyr)
library(tibble)
library(stargazer)
library(dplyr)
library(haven)
library(ggplot2)
library(zoo)

#########################################################################
# REPLICATION FILES: MILITARIZATION AND PERCEPTIONS OF LAW ENFORCEMENT  #
# JOURNAL: BJPS                                                         #
# FIGURES IN SUPP. APP. FOR:                                            #
# LEVELS OF FAVORABILITY USING MARGINAL MEANS                           #
#########################################################################

#### STEP 1: LOAD RESULTS                     ####

ratings_age <-  read.csv("MM_age.csv")
ratings_edu <-  read.csv("MM_Education.csv")
ratings_gen <-  read.csv("MM_Gender.csv")
ratings_hom <-  read.csv("MM_Homicides.csv")
ratings_ideo <-  read.csv("MM_Ideology.csv")
ratings_inc <-  read.csv("MM_Income.csv")
ratings_perc <-  read.csv("MM_Insecurity.csv")
ratings_trust <-  read.csv("MM_Trust.csv")
ratings_vic <-  read.csv("MM_Victim.csv")
ratings_pid <-  read.csv("MM_PID.csv")


#### STEP 2: ORDER FACETS                     ####

ratings_age$model<- as.factor(ratings_age$model)
ratings_age$model<- factor(ratings_age$model, levels = c("Effectiveness", "Civil Liberties",
                                                         "Corruption", "Neighborhood"))
ratings_edu$model <- factor(ratings_edu$model, levels = c("Effectiveness", "Civil Liberties",
                                                          "Corruption", "Neighborhood"))
ratings_gen$model <- factor(ratings_gen$model, levels = c("Effectiveness", "Civil Liberties",
                                                          "Corruption", "Neighborhood"))
ratings_hom$model <- factor(ratings_hom$model, levels = c("Effectiveness", "Civil Liberties",
                                                          "Corruption", "Neighborhood"))
ratings_ideo$model <- factor(ratings_ideo$model, levels = c("Effectiveness", "Civil Liberties",
                                                            "Corruption", "Neighborhood"))
ratings_inc$model <- factor(ratings_inc$model, levels = c("Effectiveness", "Civil Liberties",
                                                          "Corruption", "Neighborhood"))
ratings_perc$model <- factor(ratings_perc$model, levels = c("Effectiveness", "Civil Liberties",
                                                            "Corruption", "Neighborhood"))
ratings_trust$model <- factor(ratings_trust$model, levels = c("Effectiveness", "Civil Liberties",
                                                              "Corruption", "Neighborhood"))
ratings_vic$model <- factor(ratings_vic$model, levels = c("Effectiveness", "Civil Liberties",
                                                          "Corruption", "Neighborhood"))
ratings_pid$model <- factor(ratings_pid$model, levels = c("Effectiveness", "Civil Liberties",
                                                          "Corruption", "Neighborhood"))

#### STEP 3: SET THEME AND LABS                    ####


theme_bw1 <- function(base_size = 13, base_family = "") {
  theme_grey(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.text.x =       element_text(size = 8, colour = "black",  hjust = .5 , vjust=1),
      axis.text.y =       element_text(size = 11 , colour = "black", hjust = 0 , vjust=.5 ), # changes position of X axis text
      axis.ticks =        element_line(colour = "grey50"),
      axis.title.y =      element_text(size = base_size,angle=90,vjust=.01,hjust=.1),
      axis.title.x =      element_text(size =11),
      legend.position="bottom",
      panel.background = element_rect(fill = "#F5F5F5"),
      strip.background = element_rect(colour="black", fill="white"),
      plot.caption = element_text(size=8)
    )
}

yylab1  <- c("Estimated Marginal Means")

model<-c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood")
meanrate <-c(0.59,0.55,0.53,0.62)
hlinemeans<-data.frame(model,meanrate)

#### STEP 4: FIGURES                     ####


#--------------- FIGURE A.13 AGE ---------------------- #

mm_byage = ggplot(ratings_age, aes(y=pe,x=var, colour=measure)) +
  geom_hline(data=hlinemeans, aes(yintercept=meanrate), colour="darkgrey",linetype = "dashed") +
  coord_flip(ylim = c(.48, .68)) +
  geom_point(size=1) +
  geom_linerange(aes(ymin=lower95,ymax=upper95,colour = measure),position="dodge", size=1) +
  scale_y_continuous(name=yylab1) +
  xlab("") +
  facet_grid(.~ factor(model, levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))) +
  scale_color_manual(name="Age:", values=c("#999999", "black"), labels=c("Over 40 years old", "18-40 years old")) +
  scale_x_continuous(breaks = 15:1,
                     labels = c("Uniform:", "Police", "Military","","Weapon:","None", "Assault Rifle", "",  "Gender:", "Female", "Male", "", "Skin Color:", "Dark", "Light"))

plot(mm_byage)

#pdf(paste("Age.pdf",sep=""),width=10,height=6) 
#mm_byage =  mm_byage + theme_bw1()
#print(mm_byage)
#dev.off()

#ggsave("Age.png", plot=mm_byage, width=10,height=6, units = "in", dpi = 300)

# ---------------- FIGURE A.14 GENDER ---------------- #

mm_bygen = ggplot(ratings_gen, aes(y=pe,x=var, colour=measure)) +
  geom_hline(data=hlinemeans, aes(yintercept=meanrate), colour="darkgrey",linetype = "dashed") +
  coord_flip(ylim = c(.48, .68)) +
  geom_point(size=1) +
  geom_linerange(aes(ymin=lower95,ymax=upper95,colour = measure),position="dodge", size=1) +
  scale_y_continuous(name=yylab1) +
  xlab("") +
  facet_grid(.~ factor(model, levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))) +
  scale_color_manual(name="Gender:", values=c("#999999", "black")) +
  scale_x_continuous(breaks = 15:1,
                     labels = c("Uniform:", "Police", "Military","","Weapon:","None", "Assault Rifle", "", "Gender:", "Female", "Male", "", "Skin Color:", "Dark", "Light"))

plot(mm_bygen)

#pdf(paste("Gender.pdf",sep=""),width=10,height=6) 
#mm_bygen =  mm_bygen + theme_bw1()
#print(mm_bygen)
#dev.off()

#ggsave("Gender.png", plot=mm_bygen, width=10,height=6, units = "in", dpi = 300)

#---------------- FIGURE A.15 EDUCATION ----------------#

mm_byedu = ggplot(ratings_edu, aes(y=pe,x=var, colour=measure)) +
  geom_hline(data=hlinemeans, aes(yintercept=meanrate), colour="darkgrey",linetype = "dashed") +
  coord_flip(ylim = c(.48, .68)) +
  geom_point(size=1) +
  geom_linerange(aes(ymin=lower95,ymax=upper95,colour = measure),position="dodge", size=1) +
  scale_y_continuous(name=yylab1) +
  xlab("") +
  facet_grid(.~ factor(model, levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))) +
  scale_color_manual(name="Education:", values=c("#999999", "black"), labels=c("High school diploma or above", "Less than high school")) +
  scale_x_continuous(breaks = 15:1,
                     labels = c("Uniform:", "Police", "Military","","Weapon:","None", "Assault Rifle", "",  "Gender:", "Female", "Male", "", "Skin Color:", "Dark", "Light"))

plot(mm_byedu)

#pdf(paste("Education.pdf",sep=""),width=10,height=6) 
#mm_byedu =  mm_byedu + theme_bw1()
#print(mm_byedu)
#dev.off()

#ggsave("Education.png", plot=mm_byedu, width=10,height=6, units = "in", dpi = 300)

# ----------------FIGURE A.16 INCOME ---------------- #

mm_byinc = ggplot(ratings_inc, aes(y=pe,x=var, colour=measure)) +
  geom_hline(data=hlinemeans, aes(yintercept=meanrate), colour="darkgrey",linetype = "dashed") +
  coord_flip(ylim = c(.48, .68)) +
  geom_point(size=1) +
  geom_linerange(aes(ymin=lower95,ymax=upper95,colour = measure),position="dodge", size=1) +
  scale_y_continuous(name=yylab1) +
  xlab("") +
  facet_grid(.~ factor(model, levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))) +
  scale_color_manual(name="Income:", values=c("#999999", "black")) +
  scale_x_continuous(breaks = 15:1,
                     labels = c("Uniform:", "Police", "Military","","Weapon:","None", "Assault Rifle", "", "Gender:", "Female", "Male", "", "Skin Color:", "Dark", "Light"))

plot(mm_byinc)

#pdf(paste("Income.pdf",sep=""),width=10,height=6) 
#mm_byinc =  mm_byinc + theme_bw1()
#print(mm_byinc)
#dev.off()

#ggsave("Income.png", plot=mm_byinc, width=10,height=6, units = "in", dpi = 300)

# ---------------- FIGURE A.17 TRUST ---------------- #

mm_bytrust = ggplot(ratings_trust, aes(y=pe,x=var, colour=measure)) +
  geom_hline(data=hlinemeans, aes(yintercept=meanrate), colour="darkgrey",linetype = "dashed") +
  coord_flip(ylim = c(.45, .72)) +
  geom_point(size=1) +
  geom_linerange(aes(ymin=lower95,ymax=upper95,colour = measure),position="dodge", size=1) +
  scale_y_continuous(name=yylab1) +
  xlab("") +
  facet_grid(.~ factor(model, levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))) +
  scale_color_manual(name="Trust in the military:", values=c("#999999", "black"), label=c("5-7","1-4")) +
  scale_x_continuous(breaks = 15:1,
                     labels = c("Uniform:", "Police", "Military","","Weapon:","None", "Assault Rifle", "", "Gender:", "Female", "Male", "", "Skin Color:", "Dark", "Light"))

plot(mm_bytrust)

#pdf(paste("Trust.pdf",sep=""),width=10,height=6) 
#mm_bytrust =  mm_bytrust + theme_bw1()
#print(mm_bytrust)
#dev.off()

#ggsave("Trust.png", plot=mm_bytrust, width=10,height=6, units = "in", dpi = 300)


#---------------- FIGURE A.18 IDEOLOGY ---------------- #

mm_byideo = ggplot(ratings_ideo, aes(y=pe,x=var, colour=measure)) +
  geom_hline(data=hlinemeans, aes(yintercept=meanrate), colour="darkgrey",linetype = "dashed") +
  coord_flip(ylim = c(.48, .72)) +
  geom_point(size=1) +
  geom_linerange(aes(ymin=lower95,ymax=upper95,colour = measure),position="dodge", size=1) +
  scale_y_continuous(name=yylab1) +
  xlab("") +
  facet_grid(.~ factor(model, levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))) +
  scale_color_manual(name="Ideology:", values=c("#999999", "black")) +
  scale_x_continuous(breaks = 15:1,
                     labels = c("Uniform:", "Police", "Military","","Weapon:","None", "Assault Rifle", "", "Gender:", "Female", "Male", "", "Skin Color:", "Dark", "Light"))

plot(mm_byideo)

#pdf(paste("Ideology.pdf",sep=""),width=10,height=6) 
#mm_byideo =  mm_byideo + theme_bw1()
#print(mm_byideo)
#dev.off()

#ggsave("Ideology.png", plot=mm_byideo, width=10,height=6, units = "in", dpi = 300)

# ---------------- FIGURE A.19 VICTIMS ---------------- #

mm_byvic = ggplot(ratings_vic, aes(y=pe,x=var, colour=measure)) +
  geom_hline(data=hlinemeans, aes(yintercept=meanrate), colour="darkgrey",linetype = "dashed") +
  coord_flip(ylim = c(.45, .72)) +
  geom_point(size=1) +
  geom_linerange(aes(ymin=lower95,ymax=upper95,colour = measure),position="dodge", size=1) +
  scale_y_continuous(name=yylab1) +
  xlab("") +
  facet_grid(.~ factor(model, levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))) +
  scale_color_manual(name="Victimization:", values=c("#999999", "black")) +
  scale_x_continuous(breaks = 15:1,
                     labels = c("Uniform:", "Police", "Military","","Weapon:","None", "Assault Rifle", "", "Gender:", "Female", "Male", "", "Skin Color:", "Dark", "Light"))

plot(mm_byvic)

#pdf(paste("Victims.pdf",sep=""),width=10,height=6) 
#mm_byvic =  mm_byvic + theme_bw1()
#print(mm_byvic)
#dev.off()

#ggsave("Victims.png", plot=mm_byvic, width=10,height=6, units = "in", dpi = 300)


# ---------------- FIGURE A.21 PERCEIVED INSECURITY ---------------- #


mm_byperc = ggplot(ratings_perc, aes(y=pe,x=var, colour=measure)) +
  geom_hline(data=hlinemeans, aes(yintercept=meanrate), colour="darkgrey",linetype = "dashed") +
  coord_flip(ylim = c(.48, .68)) +
  geom_point(size=1) +
  geom_linerange(aes(ymin=lower95,ymax=upper95,colour = measure),position="dodge", size=1) +
  scale_y_continuous(name=yylab1) +
  xlab("") +
  facet_grid(.~ factor(model, levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))) +
  scale_color_manual(name="Perceptions of insecurity:", values=c("#999999", "black")) +
  scale_x_continuous(breaks = 15:1,
                     labels = c("Uniform:", "Police", "Military","","Weapon:","None", "Assault Rifle", "", "Gender:", "Female", "Male", "", "Skin Color:", "Dark", "Light"))

plot(mm_byperc)

#pdf(paste("Perception.pdf",sep=""),width=10,height=6) 
#mm_byperc =  mm_byperc + theme_bw1()
#print(mm_byperc)
#dev.off()

#ggsave("Perception.png", plot=mm_byperc, width=10,height=6, units = "in", dpi = 300)


# ---------------- FIGURE A.22 HOMICIDES ---------------- #

mm_byhom = ggplot(ratings_hom, aes(y=pe,x=var, colour=measure)) +
  geom_hline(data=hlinemeans, aes(yintercept=meanrate), colour="darkgrey",linetype = "dashed") +
  coord_flip(ylim = c(.48, .68)) +
  geom_point(size=1) +
  geom_linerange(aes(ymin=lower95,ymax=upper95,colour = measure),position="dodge", size=1) +
  scale_y_continuous(name=yylab1) +
  xlab("") +
  facet_grid(.~ factor(model, levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))) +
  scale_color_manual(name="Homicide Rate:", values=c("#999999", "black"), labels=c("15.8 or above", "Below 15.8")) +
  scale_x_continuous(breaks = 15:1,
                     labels = c("Uniform:", "Police", "Military","","Weapon:","None", "Assault Rifle ", "", "Gender:", "Female", "Male", "", "Skin Color:", "Dark", "Light"))

plot(mm_byhom)

#pdf(paste("Homicides.pdf",sep=""),width=10,height=6) 
#mm_byhom =  mm_byhom + theme_bw1()
#print(mm_byhom)
#dev.off()

#ggsave("Homicides.png", plot=mm_byhom, width=10,height=6, units = "in", dpi = 300)


# ---------------- FIGURE A.20 PID ---------------- #

# Figures SA. X

mm_bypid = ggplot(ratings_pid, aes(y=pe,x=var, colour=measure)) +
  geom_hline(data=hlinemeans, aes(yintercept=meanrate), colour="darkgrey",linetype = "dashed") +
  coord_flip(ylim = c(.43, .75)) +
  geom_point(size=1) +
  geom_linerange(aes(ymin=lower95,ymax=upper95,colour = measure),position="dodge", size=1) +
  scale_y_continuous(name=yylab1) +
  xlab("") +
  facet_grid(.~ factor(model, levels=c("Effectiveness", "Civil Liberties", "Corruption", "Neighborhood"))) +
  scale_color_manual(name="Party ID:", values=c("#C0C0C0", "#808080", "black")) +
  scale_x_continuous(breaks = 15:1,
                     labels = c("Uniform:", "Police", "Military","","Weapon:","None", "Assault Rifle", "", "Gender:", "Female", "Male", "", "Skin Color:", "Dark", "Light"))

plot(mm_bypid)

#pdf(paste("PID.pdf",sep=""),width=10,height=6) 
#mm_bypid =  mm_bypid + theme_bw1()
#print(mm_bypid)
#dev.off()

#ggsave("PID.png", plot=mm_bypid, width=10,height=6, units = "in", dpi = 300)


